#include <algorithm>
#include <string>
#include <sstream>
#include <iostream>
#include <vector>
#include <set>
#include <unistd.h>
#include <math.h>
#include "metrics.cpp"

const int MONTE_NUM = 10000;
std::vector<float> MONTE_VECTOR (MONTE_NUM, 0.0);

void generate_random_mi(int length, int n1, int n2) {
    std::set<int> s1;
    std::set<int> s2;
    for (int i=0; i<n2; i++) {
        s2.insert(i);
    }
//    std::set<int> s2 (v2.begin(), v2.end());

    srand ( unsigned ( time (NULL) ) ); 
    for (int i=0; i< MONTE_NUM; i++) {

        while (s1.size() < n1) {
            s1.insert(rand() % length);
        }

        std::vector<int> gene_ids_set;
        std::insert_iterator<std::vector<int> > gene_ids_set_ins(gene_ids_set, gene_ids_set.begin());
        std::set_intersection(s1.begin(), s1.end(), s2.begin(), s2.end(), gene_ids_set_ins);
        int n12 = gene_ids_set.size();
        MONTE_VECTOR[i] = mutual_information(length, n1, n2, n12);
        s1.clear();
    }
    std::sort(MONTE_VECTOR.begin(), MONTE_VECTOR.end());
}

float get_quantile(float level) {
    return MONTE_VECTOR[int(MONTE_NUM*level)-1];
}


int main( int argc, char *argv[] ) {

    int genomes_num = 0;
    int filter_int = -1;
    int nodes_num = -1;

    // parse options
    const char *optString = "g:f:n:";
    int opt = 0;
    opt = getopt( argc, argv, optString );
    while( opt != -1 ) {
        switch( opt ) {
            case 'g':
                genomes_num = atoi(optarg);
                break;
            case 'f':
                filter_int = atoi(optarg);
                break;
            case 'n':
                nodes_num = atoi(optarg);
                break;
            default:
		break;
        }        
        opt = getopt( argc, argv, optString );
    }
    
    std::cout << "# n1\tn2\t0.80\t090\t0.95\t0.97\t0.99" << std::endl;
    for (int n1=1; n1<=(genomes_num+1)/2; n1++) {
	for (int n2=n1; n2<=(genomes_num+1)/2; n2++) {
	    if ((nodes_num < 0) or (filter_int < 0) or (n1 % nodes_num == filter_int)){
		generate_random_mi(genomes_num, n1, n2);
		std::cout << n1 << "\t" << n2 << "\t" << \
		get_quantile(0.80) << "\t" << \
		get_quantile(0.90) << "\t" << \
		get_quantile(0.95) << "\t" << \
		get_quantile(0.97) << "\t" << \
		get_quantile(0.99) << \
		std::endl;
	    }
	}
    }
    
    return 0;
}
